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ABSTRACT 


A  first-order  stochastic  difference  equation  with  random 
coefficients  is  shown  to  have  a  solution  which  makes  the  marginal  distribution 
of  the  stationary  sequence  generated  by  the  equation  a  convex  mixture  of 
two  exponential  distributions.  This  Markovian  process  should  be  broadly 
applicable  in  stochastic  modelling  in  operations  analysis.  Moreover  it 
is  extended  quite  simply  to  a  mixed  exponential  process  with  mixed  pth-order 
autoregressive  and  qth-order  moving  average  correlation  structure.  Coupling 
of  the  processes  to  model  multivariate  situations  is  also  discussed. 
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1.  Introduction 

Sequences  of  independent,  exponentially  distributed  random  variables 
{E^}  play  a  central  role  in  stochastic  modelling  in  operations  analysis, 
either  to  characterize  the  times-between-events  in  homogeneous  Poisson  arrival 
processes  or  to  characterize,  for  example,  time  series  of  successive  service 
times  in  a  queue  or  successive  response  times  at  a  computer  terminal.  Un¬ 
fortunately  both  the  assumption  of  independence  and/or  the  assumption  of  an 
exponential  distribution  can  be  tenuous  in  practice.  The  distributional 
assumption  has  long  been  relaxed,  as  in  the  renewal  model  for  arrival  processes, 
but  independence  has  been  found  difficult  to  relax  in  simple  tractable  models. 
For  instance,  the  Box-Jenkins  or  ARMA  models, well  known  in  time  series  analysis, 
are  mainly  suitable  for  modelling  time  series  with  marginal  Gaussian  distri¬ 
butions  and  are  inappropriate  for  positive  variables,  such  as  response  times. 
However,  a  tractable  and  flexible  time  series  model  has  recently  been  intro¬ 
duced  by  Lawrance  (1980),  Lawrance  and  Lewis  (1981),  but  developed  only  for 
time  series  with  exponential  marginal  distributions.  The  development  here  is 
to  extend  this  new  model  into  mixed  exponential  variables.  Mixed  exponential 
variables  are  overdispersed  relative  to  exponential  variables  while  retaining 
their  nonmodal  distributional  aspect.  Further,  they  are  often  used  when  data 
suggests  that  exponential  assumptions  are  incorrect. 


2.  The  Model  to  be  Developed 

A  time  series  of  equally  spaced  observations  will  be  represented 
by  the  random  variables  {X^,  n  =  0,  +1,  +2,...}.  We  shall  be  concerned 
with  the  following  general  model. 


X  *  e  + 
n  n 


8  X^_  ^  w . p .  a 

0  w.p.  1-a 


(2.1) 


e  +  BV  X  .  , 

n  n  n-1  ’ 


(2.2) 


where  the  V  are  i.i.d.  binary  variables  with  P(V  ■  1)  =  1  -  P(V  =0)  =  a, 
n  J  n  n 

and  {c^}  are  assumed  to  be  an  i.i.d.  sequence.  This  was  the  two-parameter 
model  developed  in  Lawrance  and  Lewis  (1981)  for  stationary  Xn's  with 
exponential  marginal  distributions,  and  there  called  the  NEAR(l)  model  (new 
exponential  autoregressive  of  first  order) ;  here  we  develop  the  special 
one-parameter  model,  obtained  by  setting  8=1  in  (2.1), 


X  =  e  +VX  , 
n  n  n  n-1 


(2.3) 


for  {Xn>  having  a  mixed  exponential  marginal  distribution.  The  first 

and  crucial  step  is  to  determine  the  distribution  of  e  in  terms  of  the 

n 

marginal  distribution  of  -fX^} ;  the  model  (2.1)  or  (2.3)  provides  no 
guarantee  that  this  can  be  achieved.  The  model  (2.3)  turns  out  to  be  very 

j* 

tractable,  while  still  preserving  a  geometric  a  correlation  structure; 
however,  it  does  restrict  sample  path  behavior  to  a  preference  for  runs 
of  rising  values.  Developments  of  the  model  in  Section  5,  6  and  7  to 
both  higher  order  autoregressive  and  higher  order  moving  average  dependency 
overcome  this  effect  to  a  considerable  degree;  the  more  general  model  (2.2)  with 
mixed  exponential  marginals  has  been  found  to  be  very  limited  by  tractability . 


The  marginal  distribution  function  for  the  mixed  exponential  vari¬ 
ables  {X }  will  be  denoted  by  Fy(x)  where 

**  A 

-x/y  -x/y 

1  -  Fx(x)  =  T^e  +  7r2e  (2.4) 

and  ir1,  ^  >  0,  +  ir2  =  1  and  y^  y2  >  0,  ux  J  u2 i  since  i^,  v0  >  0  and 

P 1  4  V2  we  are  considering  only  a  convex  mixed  exponential  and  not  including  the 
ordinary  exponential. 

Equations  (2.3)  and  (2.4)  constitute  a  model  which  will  be  called 
NMEAR(l),  synonymous  for  'new  mixed-exponential  first-order  autoregressive 
process".  The  necessary  choice  for  the  distribution  of  is  given  in  Section  3. 
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A  mixed  exponential  random  variable  with  distribution  function 


given  by  (2.4)  will  be  denoted  by  ;  it  can  be  usefully 

written  as 


X  =  K  E' 
n  n  n 


(2.5) 


where  P(K  =  p  )  =  1  -  P(K  =  p  )  =  it,  and 
n  1  n  2  1 

random  variable.  The  customary  measure  of  its  dispersion  relative  to  an 


En  is  a  unit  exponential 


exponential  variable  is  the  coefficient  of  variation  C(X)  and 


C(X)  =  sd 


[*]- 


22  2  1/2 
sd(X)  _  {Vl  +  V2  +  7TlTr2(yrtJ2)  } 


E(X) 


(iTjP  ^  +  y21T2^ 


(2.6) 


This  is  always  greater  than  one,  its  value  in  the  exponential  case. 


3.  The  Innovation  Process  for  the  NMEAR(l)  Process 

To  establish  the  existence  and  distributional  form  of  the  innovation 
sequence  in  (2.3)  we  prove  the  following  result: 


Theorem.  For  0  _<  a  <  1  the  stationary  Markovian  sequence  defined  by 

(2.3)  has  a  convex  mixed  exponential  marginal  distribution  (2.4)  if  and  only 

if  the  i.i.d.  sequence  {e^}  has  the  convex  mixed  exponential  distribution 

function  F  (x)  given  by 
e 


-x/Yi 


-x/y. 


1  -  F£(x)  =  n^e 


+  n2e 


(x  >  0  )  , 


(3.1) 


where  r^,  n2  >  0,  n1  +  n2  =  1,  Y y2  >  0  and 


(3.2) 

Yl’  Y2  “  2  *V1  +  ^2  "  ap  -  +  y2  -ctv^2  “  4(l-a)p1P2]}  ;  (3.3) 

(3.4) 

(3.5) 


V  =  ^1^1  +  7T2v2  ' 


Y0  =  Vl  +  1,lP2  5 


ni "  (YrV/(YrY2);  n2 (y2-y0)/(y2-y1)  . 

To  prove  the  Theorem  we  need  the  following  preliminary  result: 
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Lemma.  For  0  _<  a  <  1  and  0  <  <  1. 

2 

(u  +  n2  -  ap)*  -  4(l-a)M1P2  >  (v1  +  M2  -  (2-a)y}  >  0  •  (3.6) 

Ue  have 

9  2  2  2 

{U1  +  p2-(2-a)v}/:  *  +  y2  ~  av0  “  4(l-a)(y1  +  P2  "  ay^v  +  v 

and  substituting  this  into  the  left-hand  side  of  (3.6)  and  cancelling  like 
terms,  the  inequality  (3.6)  is  seen  to  be  equivalent  to 

2 

y(Vl  +  M2)  >  PjV2  +  p  or  p(p.^  +  u2  -  p)  >  p^. 

Using  (3.2)  for  p  we  get 

[tt1p1  +  (1-^1)p2J[p1  +  P2  -  <7r1n1  +  (1-tt1)p2)]  >  PXP2  , 

[lTlWl  +  (1-,r1>»j2^(1-7ri)pi  +  nlV2^  >  U1U2  ’ 

[(Pi  -  P2)1Ti  +  }J2^^P2~P1^7T1  +  yl^  >  ^1^2  * 
or 

2  2  2 

-  (p1~p2)  Wj.  +  (p1-y2^  ’’i  +  ^1^2  >  ^1^2 
or 

0  <  7T1  <  1 

which  is  true  by  assumption.  Thus  the  Lemma  is  proved. 

Proof  of  the  Theorem 

For  the  stationary  sequence  given  by  (2.3)  we  define 

$x(t)  -  E{ exp(-tX)}  ,  <f> e (t)  =  E{ exp(-te) }  ,  (3.7) 

so  that  transforming  both  sides  of  (2.3)  gives 

<J>£(t)  ■  4>x(t)/{ai<ix(t)  +  (1-a) }  .  (3.8) 

With  $x(t)  *  (1  +  Y0t)/{(1  +  P1t)(l  +  p2t)}  from  (2.4),  (3.8)  yields 
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$e(t)  “  (1  +  YqO/U  +  (v1  +  V2  -  ay)t  +  (l-cOy^t2}  .  (3.9) 


The  reciprocal  roots  of  the  quadratic  denominator  are  given  by  (3.3)  and  are 
real  by  the  Lemma;  furthermore  they  are  distinct  and  nonnegative  since 
4(l-a)vjP2  >0.  A  partial  fraction  expansion  of  (3.9)  then  gives 

4»e  ( t)  =  n1(l  +  y^t)  1  +  n2(l  +  Y2t)  1  •  (3.10) 

This  inverts  to  the  required  mixed  exponential  distribution  provided  that 
the  weights  and  n2  are  both  positive  and  sum  to  one;  the  latter  is 

clearly  true  from  (3.5).  Now  >  0  requires,  from  (3.5),  that  y^  >  Yq> 
and  n2  >  0  requires  that  y2  <  Yq-  The  first  of  these  is  equivalent 
to 

\  {Vl  +  u2  -  ay  +  /[(u1  +  P2  -  avi)2  -  4(1-0^^]}  >  v1  +  v2  ~  P 
or 

/[(yx  +  y2  -  ay)2  -  4(1-0^^]  >  y]L  +  y2  -  (2-a)y  (3.11) 

and  the  second  to 

2 

/{(Pj^  +  y2  -  ay)  -  4(l-a)y;]y2}  >  -  {y1  +  V2  ~  •  (3.12) 

Squaring  (3.11)  and  (3.12)  show  that  the  condition  for  Yj_  >  Yq  and 
Y2  <  Yq  is  just  the  condition  proved  in  the  Lemma.  The  key  requirement 
is  that  should  be  a  non-zero  probability. 


4.  Utility  of  the  Theorem 

(i)  Note  that  the  required  distribution  for  the  i.i.d.  e  's  is 

n 

again  the  convex  mixed  exponential;  the  e^'s  can  then  be  generated  as 


n 


where  P(L 


V 


1  -  P(L 


V 


L  E 
n  n 


n,  and  -  E 
1  n 


(4.1) 

is  a  unit-mean  exponential. 


n  1  n 

Thus  the  model  (2.3)  gives  a  means  of  probabilistically  transforming  an  i.i.d. 

exponential  sequence  into  a  first  order  autoregressive  Markov  mixed  exponential 
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sequence  with  autocorrelations  p  -  ar  >  0.  If  a  =  0.  {X  }  are  i.i.d. 

r  —  n 

ME(it1,y1;ir2,y2)  random  variables.  Since  the  process  is  Markovian  it  is  easy 

to  show  that  setting  Xq  to  be  a  ME(Tr^,y^;ir2,p2)  random  variable  independent 

°f  el*e2* ’ *  *  makes  the  Process  {Xn,  n  =  1,2,...}  stationary.  Finally  it 

is  easy  to  see  that  the  regression  of  X  on  X  ,  is  linear:  from  (2.3) 

n  n-1 

E(X  |X  .«x)  -  E(e  )  +  ax. 
n  n-1  n 


(ii)  The  theorem  also  holds  when  the  marginal  distribution  of  X 


is  the  zero-jump  exponential.  In  our  notation,  let  y^  =  0  and  y2  =  y', 
so  there  is  a  probability  that  X  *  0.  This  structure  is  related  to  the 

queuing  situation  where  the  waiting  time  can  be  either  zero  when  no  customer 
is  being  served,  or  an  exponential  length  of  time  when  the  server  is  already 
busy.  The  theorem  gives  the  following  form  of  e  ; 


I 


(1-aOy  'E 
z  n 


w.p.  (l-a)7r2/(l-air2) 


w.p.  if1/(l-c«r2) 


(4.2) 


Thus  the  e  variable  also  has  a  zero-jump  exponential  distribution, 
n 

(ill)  There  are  two  cases  in  which  a  mixed  exponential  marginal 
distribution  ME^.y^j^.y^  of  {Xn>  can  reduce  to  a  single  exponential 
distribution;  these  are  when  y^  =  y2  and  when  either  tt^  or  ir2  is  zero. 
Neither  case  is  covered  by  the  theorem,  but  each  can  be  treated  directly 
from  the  model  (2.3).  With  y  as  the  mean  of  the  exponential  { 
sequence,  an  exponential  sequence  with  mean  (l-a)y.  This 

situation  is  the  TEAR(l)  process  studied  in  Lavrance  and  Lewis  (1981). 


5.  Mixed  Exponential  Moving  Average  Processes 

It  is  possible  to  think  of  (2.3)  and  its  mixed  exponential  solution 
as  a  way  of  combining  independent  ME(n^,Y^;n2 ^2)  and  ME(iT^,y^;7r2,u2) 

variables  into  another  ME(tt^,u  ^;7r2  ,y2)  variable.  This  key  result  allows 
the  constructing  (cf.  Lawrance  and  Lewis,  1977)  of  a  first  order  moving 


average  process 
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n 


L  E  +  V  K  ,E  , 
n  n  n  n-1  n-1 


0,  +  1,  +  2, 


(5.1) 


where  V^,  and  are  defined  at  (2.2),  (2.5)  and  (4.1),  respectively. 

The  dependency  parameter  Is  a.  By  the  results  of  Section  3,  X^  will  have 

the  ME(ir^, tt2»  U2^  distribution  (2.4),  but  the  dependency  will  be  that  of 

a  first  order  moving  average  process;  the  model  will  be  denoted  as  NMEMA(l) . 

The  moving  average  aspect  is  clear  because  Xn  and  xn+^  have  just  E^ 

in  common  while  X  and  X  (r  >  1)  will  have  no  E  's  in  common  and  thus 
n  n+r  n 

be  independent.  However,  the  NMEMA(l)  process  has  a  very  limited  range  of 

correlations;  corr(X  ,X  ,,)  is  easily  derived  as 

n  n+l 

2  2  2  2 

P-L  =  a(l-a)(TT1y1  +  Tr2y2)  ZU-^  +  ir2P2  +ir1ir2^i“y2^  ^  (5.2) 

=  a(1-'t/C2(X)  .  (5.3) 

The  a(l-a)  is  always  less  than  one-fourth  and  the  coefficient  of  variation 
C(X)  of  the  mixed  exponential  distribution  is  always  greater  than  one;  thus 
an  overall  bound  below  one-fourth  is  implied .  Note  that  the  maximum  allowable 
value  of  decreases  as  C(X)  increases,  i.e.  as  the  dispersion  of  X 

relative  to  an  exponential  random  variable  for  which  C(X)  =  1,  increases.  A 
stationary  sequence  is  clearly  obtained  for  n  =  1,2,...  by  starting  the 
moving  average  over  Eq ,  E^ . 


6.  Higher-Order  Autoregressive  Mixed  Exponential  Models 

Although  it  is  possible  to  construct  higher-order  analogs  of  the 
autoregressive  equation  (2.1),  as  in  Lawrance  and  Lewis  (1980),  a  simple 
p-th  order  extension  of  the  model  (2.3)  is  obtained  by  letting  xn  ^  be  a 
probability  mixture  of  xn_^»  ’  xn  p’  the  xn's  are  mar8*naHy 

ME(7r^,y^;ir2»y2)  then  this  mixture  will  be  also,  and  the 

error  sequence  in  the  autoregressive  extension 
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0 


w.p.  1  -  (c^  +•*•+  otp) 


E  + 

n 


rx  - 

n-1 


n-2 


w.p.  a2 
w.p.  a2 


(6.1) 


X 

-P 


w.p.  ct 

P 


Is  given  as  at  (3.1)  to  (3.5)  to  be  ME(n^,y^;n2 > Y2)  with  a  replaced  by 
a.  +  a„  +•••+  a  .  The  autocorrelations  now  follow  pth  order  linear,  constant 

1  2  p 

coefficient  difference  equations  paralleling  those  of  the  standard  auto¬ 
regressive  models. 

The  second-order  situation  NMEAR(2)  is  probably  as  high  as  one 
would  want  to  go  in  modelling  data,  in  particular  to  take  care  of  situations 
where  the  marginal  distribution  is  mixed  exponential  but  either  the  process 
is  not  first  order  Markov  or  the  sample  paths  do  not  tend  to  run 
up.  In  this  case  simple  computations  give 

P1  =  a]/ (l-c»2^  2  «1  *  p2  =  a2  +  alPl  '  (6.2) 

The  range  of  p^  and  P2  combinations  is  restricted  by  0  <  a^,  02  (  1,  and 
«1  +  a2  <  1.  In  particular  any  first  autocorrelation  p*  between  zero  and  one 
is  attainable  for  values  of  and  satisfying  =  1-a^/p*.  Since  this 

is  linear,  ranges  between  0  and  p*  as  ranges  from  1  to  0  , 

where  a2  =  0  gives  the  NMEAR(l)  case.  Note  that  ^  =  0  implies  that 
*X1’  are  independent  of  {X2,  X^,...}  and  that  the  odd  and  the  even 

sequences  are  governed  by  separate  but  identical  first  order  processes.  The 
effect,  practically,  of  being  able  to  choose  a  small  for  a  given  p 

will  be  to  adjust  run-up  behavior. 
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7.  Mixed  Autoregressive-Moving  Average  Extensions 

It  is  possible  to  combine  the  autoregressive  structure  (6.1) 
and  the  moving  average  structure  (5.1),  as  was  done  by  Jacobs  and  Lewis  (1977) 
and  Lawrance  and  Lewis  (1980)  for  exponential  variables,  to  obtain  a  complete 
NMEARMA(p.q)  process  with  mixed  pth  order  autoregressive-qth  order  moving 
average  correlation  structure.  This  yields  a  much  richer  process  in  terms 
of  sample-path  behavior  but  since  the  process  is  not  Markovian  it  is  diffi¬ 
cult  to  do  estimation  for  the  parameters.  Thus  for  the  sake  of  completeness 
we  consider  only  the  (1,1)  model,  extending  the  allowable  range  of  marginal 
behavior  of  the  EARMA(1,1)  process  of  Jacobs  and  Lewis  (1977);  higher-order 
extensions  are  obvious.  The  NMEARMA(1, 1)  process  is  defined  by  the  equations 


X  =  L  (a) E  +  V  (a)Y  .  , 
n  n  n  n  n-1 


Y  -  L'(a')E  +  V’(a')Y  , 

n  n  n  n  n— 1 


n=0,  +1,  +2,  .... 


(7.1) 


where  the  independent  i.i.d.  sequences  V  (a)  and  V'(a')  have 

n  n 

P{V  (a)  *  1}  =  a,  PtV^a')”  1}  =  a’.  Also  L  (a)  and  L'(a')  are  inde- 
n  n  n  n 

pendent  with  distributions  given  at  (4.1).  By  the  Theorem  of  Section  3  the 

(X  }  sequence  has  a  ME(ir1  ,n,  ,p_)  marginal  distribution  and  a  correlation 

n  llzz 

structure  typical  of  standard  ARMA(1,1)  models.  Other  properties  of  the 
model  go  through  in  complete  analogy  with  the  EARMA(1,1)  process.  In  fact 
they  will  be  simpler  since  that  process  is  based  on  the  EAR(l)  model  of 
Gaver  and  Lewis  (1980)  whose  zero-defect  gives  most  of  the  problems  in 
the  limit  theorems  of  Jacobs  and  Lewis  (1977). 


8.  Modelling  and  Multivariate  Extensions 

The  mixed  exponential  time  series  generated  by  (2.3)  and  written. 


using  (4.1)  and  (2.5),  as 


X  *  L  E  +  V  X  , , 
n  n  n  n  n-1 


n  =  0,  +  1,  +  2, . 
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is  a  simple,  random  linear  combination  of  random  variables.  As  with  the 
ordinary  ARMA  models,  this  makes  it  ideal  for  modelling  complex  systems  and 
for  extension  to  multivariate  time  series.  The  use  of  EARMA  models  in 


modelling  queues  has  been  discussed  in  Lewis  and  Shedler  (1977)  and  Jacobs 

(1980);  the  NEAR(l)  structure  used  here  is  actually  better  for  this  purpose 

since,  unlike  the  EAR(l)  structure,  the  error  term  does  not  disappear. 

These  applications  of  the  NMEAR(l)  process  will  be  discussed  elsewhere. 

Here  we  only  discuss  direct  multivariate  extensions. 

Thus  let  (E^,  E")  be  a  bivariate  pair  with  unit  exponential 

marginals  and  let  X'  and  X"  be 

n  n 


X’ 

n 


L'E' 
n  n 


V'X* 
n  n- 


1  * 


X”  =  L"E"  +  V"X"  . 
n  n  n  n  n-1 


n  “  °»  ±  ±  2. • • •  (8.1) 


where  L'  is  chosen  to  make  X'  be  ME(x '  ,p '  ;ir '  ,p ' )  and  L"  is  chosen  to 
n  n  1  I  z  z  n 

make  X"  be  ME(7r',,y";7r«,p^) .  The  i.i.d.  squences  (L',L">  and 

n  X  i.  4  z  n  n  n  n 

could  be  dependent  within  pairs.  The  sequences  { X^}  and  f X”)  will  have 

negative  cross-correlation  if  the  {E',E")  pairs  are  negatively  correlated. 

n  n 

A  cross-coupled  version  of  this  process,  as  in  Gaver  and  Lewis  (1980),  will 
create  alternating  and  posibly  negative  autorcorrelation  in  and  Xjj. 

Many  other  multivariate  extensions  are  possible. 


9.  Discussion 

This  paper  has  focussed  on  replacing  the  'Poisson  process'  assumptions  of 
exponentiality  and  independence  which  are  often  made  in  modelling  sequences 
of  identically  distributed  positive  random  variables.  The  exponential 
assumption  has  been  replaced  by  a  mixed  exponential  marginal  distribution, 
allowing  more  variable  behavior,  while  the  independence  has  been  replaced 


by  autogressive  and  moving  average  dependence.  Utility  of  the  model  has  been 
referred  to  in  queuing  and  multivariate  situations.  Simulation  of  the  model 
la  possible  from  easily  obtainable  i.i.d.  exponential  sequences. 

Various  aspects  of  the  model  are  currently  being  explored  or  are  open 
to  development.  Extensions  incorporating  negative  dependency  can  be  treated 
by  the  method  of  cross-coupling  used  in  Gaver  and  Lewis  (1980)  and  Lawrance 
and  Lewis  (1981).  Estimation  via  the  likelihood  is  available  in  the  pure  auto¬ 
regressive  models,  but  the  standard  assumptions  are  not  satisfied  because  of 
discontinuities  in  the  likelihood.  This  aspect  is  investigated  in  Raftery 
(1980a, b)  under  the  assumption  of  an  exponential  marginal  distribution. 
Extensions  to  nonconvex  mixed  exponential  solutions  may  be  useful  when  the 
marginal  distribution  is  less  dispensed  relative  to  the  exponential  distri¬ 
bution;  one  such  case  is  from  the  sum  of  two  independent  non- identical  expo¬ 
nentials.  However,  the  model  has  no  Gamma  solution.  A  more  flexible  two- 
parameter  autoregressive  model  is  available,  but  its  development  is  limited 
by  mathematical  tractability. 
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